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Abstract 

We consider the short-time energy relaxation of the dilute SK model. We show that the more 
modular the system, the more rapidly the energy decays at short times. Conversely, a more modular 
system reaches a less favorable energy at long times in a static environment. We use these results 
to discuss the dynamics of the modularity order parameter in a system for which the coupling 
parameters of the dilute SK model change in time, due to a changing environment. Modularity 
endows the spin glass with a better response function in a changing environment. 

PACS numbers: 87.10.-e,87.15.A-,87.23.Kg,87.23.Cc 
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INTRODUCTION 



Biological systems are modular, and the organization of their genetic material reflects this 
modularity ^^iM- Genes are organized into exons, and expression of different exons allows 
one gene to produce multiple proteins. In many cases, each exon confers a distinct function 
to the protein, and so an exon is a modular element. Similarly, while genes interact, each 
gene confers function or functions to the organism, and so a gene is a modular element. 
Collections of related genes clustered together on the genome occasionally form a functional 
unit termed an operon, which is also a module. 

Complementary to this modularity is a set of evolutionary dynamics that evolves the 
genetic material of biological systems. Horizontal gene transfer (HGT) is a mechanism by 
which genes, pieces of genes, or multiple genes are transferred from one individual to another. 
The two individuals need not be of the same species. For two individuals of the same species, 
recombination is possible, and recombination may be viewed as a subset of HGT. 

The organization of biology into modules simultaneously restricts the possibilities for 
function, because the modular organization is a subset of all possible organizations, and may 
lead to more rapid evolution, because the evolution occurs in a vastly restricted modular 
subspace of all possibilities js, 6|. There is, then, a tension between the increased ability to 
evolve that modularity confers upon a system and the penalty that modularity imposes upon 
a system. The amount of modularity that evolves in a system reflects a balance between 
these two competing effects, as a function of the timescale at which selection occurs jsj. 

The degree of modularity that evolves in a system reflects the amount of pressure upon 
the system to evolve 0]. For a system in a changing environment, the rate of growth of 
modularity from an initially non-modular state is roughly proportional to the frequency of 
environmental change sl. Similarly, if there is selection pressure for a system to evolve 
rapidly to steady state from an initially unfit state, modularity can arise |8|. A changing 
landscape essentially selects for the response function of the system, at the time scale of 
the environmental change [sj . The coupling between transport and a heterogeneous spatial 
environment also leads to the emergence of modularity 6|. 

We here analyze a simpler model than biological evolution, that of energy relaxation 
in a spin glass. The model is similar in spirit to the spin-glass models that have been 
introduced to analyze the relation between genotype and phenotype evolution 
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body contributions to the fitness function in biology, leading to a rugged fitness landscape 
and glassy evolutionary dynamics, are increasingly thought to be an important factor in 
evolution 13|. That is, biological fitness functions may be characterized as instances of 



fitness functions taken from a spin glass ensemble. Importantly, though, biological fitness 
functions have a modular structure, and their dependence on the underlying variables is 



somewhat separable 
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lq |. Glassy evolutionary dynamics has been noted a number of 



We here analyze, within the context of statistical mechanics rather than a detailed evo- 
lutionary model, the dependence of a spin glass response function on modularity of the 
interactions. We carry out these calculations for short times, to make predictions for how 
a spin glass would relax upon change of its coupling parameters. In the present context, a 
change of environment means a change of the coupling parameters in the spin glass. Nu- 
merical calculations have shown that the energy per spin relaxes at different rates for spin 
glass systems of different sizes [l^, and these simulations provide additional motivation for 
the present calculations. 

In a full replica symmetry breaking calculation, we will show that the response func- 
tion at short times depends on the modularity. Since modularity has been argued to be 
a relevant, emergent order parameter 

mm 

, we consider the ensemble of spin glass 
Hamiltonians parametrized by modularity, M. Near the replica symmetry breaking tran- 
sition, where analytic calculations are possible, we will show that the response function 
increases monotonically with modularity. This calculation contains two technical advances: 
a generalization of the dynamical equations of magnetization and energy [2^ to the dilute 
SK model and a full replica symmetry breaking calculation to determine the form that these 
equations take near the spin glass phase transition. 



MODEL 



The focus of the present study is how to introduce modularity to the SK model, and the 
resulting short-time dynamics. The coupling matrix must have local structure, and it must 
be sparse, as modularity can not be identified in a fully connected matrix. The non-zero 
entries in coupling matrix are shown in Fig. [T] 

We here consider a spin glass model that generically incorporates sparseness and mod- 
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FIG. 1: Shown is a simplified view of the couplings in the dilute SK model. In this figure, we 
consider a system of size = 20. If spin i interacts with spin j, a dot is displayed at matrix 
position Each position i interacts on average with C other positions. Here C = 6. Left) A 
non-modular structure, M = 0. Middle) A moderately modular structure, M = 2/3. Right) A 
fully modular structure, M = 1. The matrix shown here is the connection matrix, denoted by the 
symbol A. Here, there are two modules, each of size L = 10. We define modularity from the excess 
number of interactions within the two L x L block diagonals over that expected based upon the 
probability observed outside the block diagonals. This number is divided by the total number of 
interactions to give the modularity, M. 

ularity. The connection matrix for a given system a is denoted by A° with A^- = 0, 1, as 
shown in Fig. [H If we were performing a fully evolutionary calculation we might call this 
system a short protein, with a fold structure indexed by a. Each spin i is connected to C 
other spins, on average. Putting these points together, our simplified model is a dilute SK 
model: 



with Jij = Jzij where z is a quenched Gaussian with zero mean and variance 1/C. The 
number C is the average number of connections, and so in the absence of modularity 
P(Ajj) = (1 — C/N)6^-^fi + {C/N)6/^-^^. We have cTj = ±1. The spin dynamics is gov- 
erned by Glauber dynamics such that the rate to flip spin k in the sequence is given by 
WkiW}) = 1(1 - o-fctanh/3/ifc) where hk = Y.j^kJkj^kj(Tj = Jzk, with Zk = Y^j^k^kj^kjCTj. 

Here we consider the spin glass ensemble to be parametrized by modularity, such that 
there is an excess of interactions in A along the N/ ki x N/ ki block diagonals of the N x N 
connection matrix. Thus, the probability of a connection is Cq when [kii/N\ ^ [_kij/N\ 
and Ci when [kii/N\ = [A;ij'/A^J. The number of connections is C = Cq + (Ci — Coj/ki. 
The modularity is defined by M = (Ci — Cq) / {kiC). To see the spin glass phase, the system 




(1) 
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must be macroscopic, N oo. In addition, the module size must be large, so that the glass 
phase appears, and so ki must be large. We require ki — )■ oo, but ki/N — )■ 0. To calculate 
some of the coefficients, we will require C be large, although C ^ A^, N/ki. 

We define the total magnetization m = Xlili '^i fitness r = —H/{JN). We 

project the microscopic probability of a given state, -Pt(cr), onto these order parameters. 
These order parameters evolve according to [2^ 

d/TTl 

= J dzDm,r]t[z]tsjih. (3Jz — m 

dv 

-r = J dzDm,r;t[z]z tanh (3 Jz — 2r (2) 

(JjL 

where 

n M - r PtjcrMm - m{a)]6[r - r(a)]^ ^jz - z.ja)] 

We assume that D^^z] .s self-averagmg ov. the disorder, wffich numerical simulations 
out to intermediate times seem to support [20|. We will also assume equipartitioning of 
probability in the macroscopic subshell (m,r) [2^. These assumptions allow us to drop 
Ptic) and to perform the averages over the quenched random Zij and Ajj variables: 

n ri r / 5{[m - m{a)]5[r - r{a)]j^Y.k=i ^[^ ' ^k{a)] 
JV^oo \ 6{[m - m{(r)]6[r - r{(T)] 

REPLICA SYMMETRY BREAKING 

We now proceed to analytically calculate the averages required to determine the solution 
to Eq. ([2]). We define w{a) = 6{[m — m{a)]6[r — r[a)]. We use the replica expression in the 
form 

Tr^w(cr) 

TT^i,„^nw{a^)<l>{cr^)w{a^) ■ ■ ■ w((t") 
~ Tr^i...^ntf;(cri) ■ ■ ■w((T") 

Tr„i „r.w(a^)^(a^)w(a^) ■ ■ ■wia'') 




lini , 

n->o [Tr<^w(cr)]" 

lim TT^i,„,nw{a^)^{a^)w{a^) ■ ■ ■ ^(a") (5) 

n— >0 
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to write 



Dm,r:t[z\ 



1 ^ 

lim lim — V (Tr^i,„^n5[z - Zk{a^)]w{a^)w{a^) ■ ■ ■ (a")) . 



},{A,,} 



(6) 



k=l 



Using the Fourier representation of the delta function, we find 20 1 



Dm r \z] = hm hm — [ — 

k=l 



n 



Ndnia Ndra 
~27r 2tt~ 



X ( e 



(7) 



{2ij}>{Aij} 

We average the quantity in brackets over the Aij, setting A; = 1 by permutation symmetry 



to find 



n 



1 _ 9l\ + ^e-"-^u'^]-iEc^^c.af2i,<T- 
N N 



V N N 

n 



c 



n 

D I- 



1 _ ^ ] + ^e-iEc^-<.'^?^..'^" 
iV / iV 



(8) 



where c,C indicates within the N/ki x N/ki block diagonals and rf, D indicates outside these 
block diagonals, lowercase indicates i = 1 and uppercase indicates i > 1. Recognizing that 
Cq/N and Ci/N are small, so that the above expression can be written in exponential form, 
Eq. ([7]) becomes 

N 



Dm r \z] = lim lim — [ — 

k=l 



n 



Ndma Ndva 



Ea Eb] ( (exp(-i E. ----f 



where due to symmetry, A indicates the upper half of the connection matrix, B indicates 
the upper half of the N/ki x N/ki block diagonals, a indicates the row with j ^ k, and b 
indicates the row of the block diagonal with j ^ k. 

We calculate these averages using replica symmetry (RS), 1-step replica symmetry break- 
ing (1-RSB), and full replica symmetry breaking (FRSB). We introduce overlap parameters 
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for the whole matrix and for the block- diagonal part of the matrix as 



N 

i 



QaiSfsio- ) = 

B f \ ^1 a 



j, block 



Csi^) = I E ^rc^f^z^f (10) 

j, block 

The four sums inside the exponential in Eq. iQ sum to Ntlj[q{a)] + glai, q{o')] where 
^[g(a)] = l[To(r~)-l] (Co + CM) 

a<0 

+... (11) 



and 



g[a^, q{a)] = {ChTo{r) - To(f)) (Co + CM) 

+ J2 (chTfir) - Tf{r)) a^^, {C,qtp{a) + CMq^p{a)) 



a<l3 

+ J2 ShT,-{r)a^ {Coqticr) + CMgfja)) 

a 

+ J2 ShTf\r)a'^a^,aJ{Coqt0,icr) + CMqt,^{a)) 

a</3<7 

+ ... (12) 

where terms higher order in the spin overlaps have been omitted. Here T, ChT, and ShT 
are combinatorial factors: 



jiai«2---«fc _ / ira^Zjj) ■ ■ ■ia,v}i{—ira^Zij) JJ^ cosh(z 

\ w=l 

ChT^^"^'""'' (f) = ( cosh{ixZij) ta,nh{—ifaiZij) ■ ■ ■ta,nh{— if a^Zij) Y\ cosh( 



'W=l 



ShT^^°''^"'°'''{r) = ( sinh(— tanh(— irQ,j2;jj) ■ ■ ■ tanh(— ira^^jj) ]^ cosh(2r^2;jj) \ (13) 



Expanding these in p and 1/C: 



= 1 + /(2C) + E/^t/(8C^) + E 3pIpI'/(4C^') + 

UI = 1 ui w<w' 

Tf = PaPp/C - PaP^/C'ipl + pD + ^p^pp/{2C') Y.pI + ... 

w 

T^^"'^ = {3/C^)papl3pyp5 + ■ ■ ■ 

ChTo = 1 - xV(2C) + xV(8C2) + . . . 

ChT2 = p^C - 2pyC^ - 3p2xV(2C2) + . . . 

ChT^ = Sp^C^ + ... 

ShTi = {-ix)p/C - {-ix)p^/C^ + {-ixfp/{2C^) + ... 

ShTs = {-ix)3pyC^ + ... 



(14) 



We find a final expression of 



Here 



Dm r\z] = lim lim 



dx -n- NdrUa Ndra 



n 



n 

/3=1 



TV-i^oo n.-5>0 } 2tT 271 2n 

a=l 

Ndq^pdqtp Nd~q%dq^^ 



27r 



27r 



X dq^i^^^dq^p^^ ^ dq^p^^dq^f^^^ c'^^c^-^ 



n 

7,(5=1 



2n 



ki - 1 



27r 



where g{a) = g[a, q{a) — > q] and 



in, out 



'Jf^g9(o")g^in,out(o-) 



(15) 



(16) 



. (3^ o</3 cr</3<7<5 



X 



out 



. (3^ o</3 a</3<7<(5 



(17) 
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and 



a a<l3 
a 

a</3<7<<5 

1 fci — 1 

+ — InTr^e^-^'") + ^- In Tr^e^°"*^'") 

ki ki 



(18) 



In the large N limit, these integrals reduce to a saddle point calculation, and for stability 
we find ifia = ifia and = ip^. 
We find 



m 



1 fci - 1 , > 

ki ki 



_d_ 

dpo 



■Hp) 



(19) 



and the overlap parameters to be the expected multipoint averages of the spins, with q^ = 
0{q^ /ki) for large ki. We now consider the zero net magnetization case, m = 0. We initiate 
with a random distribution of spins and watch the relaxation to equilibrium. Equation ( !T9|) 
for the fitness becomes 



b2 

1/9 



3p3 



+0(pM/c^) 



E ((i-^)47/ + ^?iV) 



l</3<7<<5 
5 1 /r^2^ 



(20) 



Note that this equation contains order parameters to all orders. Near the phase transition, 
we will keep terms to third order in p. 



REPLICA ANALYSIS 



With full replica symmetry breaking, the matrix qaf} is equal to q{x), where x is the depth 
of the tree at which a and (3 are connected. The fitness equation becomes 



2r 



p-c 



il-M) / q'^ixydx 



-M [ q^{xfdx 
Jo 



+ 0{p') 



(21) 



Defining the dynamical fitness as / = —\imn^of/n, we find 



1 2 , 1 2r 



/ = -ln2 + pr-^p^ + ^p^[(l-M)g^(l) + Mg^(l)] 



— -p" / dx 







[(l-M)/'(x)+Mg^'(a;) 



lim — 

n-i>0 n 



1 

— ln(JJcoshpM«);f 



/ti - 1 

H ; 1^(11 P^a)n 



out 



(22) 



where {uaUpY" = {1 — M)q^i^ + kiMq^^ and {uaUp)'^^^ = (1 — M)q^r ■ For replica symmetry, 



the g-dependent terms in Eq. (l22l) reduce to Eq. (2.31) in 21|, for 1-step replica symmetry 
breaking they reduce to Eq. (3.30), and for FRSB near Vc evaluating the traces up to fourth 
order in q they reduce to Eq. (3.45) with /3^J^ = p^kiM. 

The critical value of p is pc = 1/ yJkiM . Near the critical value of p for which the spin 
glass phase appears, we find g'^(x) = and 



x/2, < X < 2gi 
gi, 2gi < X < 1 



with qi = gRs = (p^fciM - l)/2. We, thus, find 



(23) 



2r = p[l-Mqlil-Aqi/3)] 



(24) 



And we find 

Dor = / {dx/2'K) exp(— x^/2 — ixz) cosh(ixp) 



1 + xVM^ / g^(x)^c/x/2 





The fitness equation becomes 

dr 
'dt 



A-2r + 45[1 + (7 - l)e(r - r^)y 



(25) 



(26) 



where rc ~ 1/ (2i//ciM) and 9 is the Heavyside step function, which implies 



r{t) = A 



t-r + 



AAB + 2 



+ 4S(7 - l)r 



x{t-t,)Q{t-Q+0[t\ {t-Q^] 



(27) 
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Defining e = {p^ — pi)/ pi we find the replica symmetric and full replica symmetry breaking 
results to be the same, 

7RS = 7FRSB = 1 + M(2 - M)eV4 . (28) 
The one-step calculation gives 

71RSB = 1 + (1 - mi)M(2 - M)eV(2 - mif , (29) 

and the multiple-step result converges to the FRSB result in the limit. The constants A and 
B are given by 

A = j dz/V2^exp{-z^/2)ztsi-n.hf3Jz[l + {z^ -6z^ + 3)/{8C) 

B = (1/2) J dz/V2^exp{-zy2)ztaBh(3Jz{z^ - 1)[1 + (^^ - 3)V(8C)] . (30) 

For a typical case of a dozen connections, C = 12, we find A = 0.597693, B = 0.481125 for 
/3J = 1 and A = 0.789573,5 = 0.452965 for /3J = oo. We, thus, find that iAB > 1 for 
T < Tc. Shown in Fig. |2] are the RS and FRSB solutions to Eq. ( l26|) . There is a kink in 
the FRSB solution at tc, and after tc the FRSB solution relaxes faster than does the RS 
solution. 



DISCUSSION 

A modular spin glass relaxes to a higher fitness value than does a non-modular spin 
glass at short times, beyond a critical time tc- The greater the modularity, the greater the 
response function. From Eq. fl27j) we see that the rate of fitness increase is larger for larger 
M. Thus, selection for systems with large short-time response function will identify modular 
systems. Another perspective is to set M = 1 so that all the connections are within the 
L X L, L = N/ki, block diagonals and see how the response depends on the block size. The 
parameter L is a measure of the effective modularity in the system, with smaller L indicating 
greater effective modularity. Since p^ = L/{NM), the spin glass transition occurs earlier 
and is larger for the system with smaller L. Thus, the system with smaller L has a larger 
fitness at short times. 

There is an optimal modularity at a given timescale. At short times, the system with 
greatest fitness has small L. At long times, the system with the greatest fitness has large 
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FIG. 2: Shown is the RS (solid) and FRSB (dashed) solution to Eq. ([26]). After the critical point, 
the FRSB solution relaxes faster than does the RS solution. Here M = 1, e = 1, and tc = 1/2; the 
theory is asymptotic in the limit of small e and tc- 



L, because more of the phase space is accessible to the connection matrix. At intermediate 
times, the optimal system will have an intermediate L, with the optimal L monotonically 
increasing with time. A system with smaller L has a less favorable infinite time fitness. 



according to r* 



22|. Arguing that the barriers to equilibration of a larger 



system further down in the Parisi hierarchy cannot be greater than the distance of the 
smaller system from the ground state, we would expect tERG ~ toexp(cL^/^) 23|-|25|. We 
expect logarithmic convergence at long time 26|. Smoothing the short time behavior, we 
expect the fitness to follow r/^(t) = r — a,L~'^^^ tanht — b[l + \n{l -\- 1 / t^RG)]''^^'^ where u = 1 
to have the expected L dependence at large time. Figure [3] shows the crossing behavior and 
illustrates the optimal systems size as a function of time. Numerical simulations exhibit the 
energy relaxation as a function of time and system size that is shown in this figure 19 1. 
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FIG. 3: Shown is the fitness an equihbrating spin glass of varying size, L (dotted=5'^ short 
dashed=6^, long dashed=7^, and solid=8'^), as a function of time. The system size with the 
highest fitness is a monotonically increasing function of time. 

SUMMARY 

We have performed a full replica-symmetry breaking calculation for the dynamics of a 
dilute SK model. Correlations in this model were defined by a connection matrix, which was 
parametrized by its modularity. We showed that energy relaxation of the dilute SK model 
is quicker for a modular connection matrix than for a non-modular connection matrix. We 
suggested that if the dilute SK model is interpreted as a rough model for evolution of 
biological structure, the present results illustrate a selective pressure for modularity to arise 
in biological populations evolving in changing environments. This interpretation rationalizes 
a number of empirical observations for increased modularity in changing environments 

Interestingly, in biology horizontal gene transfer significantly enhances the emergence 
of modularity jsj. In statistical mechanics, the present calculation shows that modularity 
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increases the short-time response function in a single physical replica. Calculation of the 
effect of horizontal gene transfer on the dynamics of modularity is an open problem. 
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